#ifndef PCPS_LOCAL_CORRELATION_HEADER
#define PCPS_LOCAL_CORRELATION_HEADER

#include <src/pcps.h>
#include <src/operator.h>
#include <src/coeff_value.h>

namespace pcps {

  //-------------------------------------------------------------------------------
  // pcps::local_correlation_functions -- Computes the contribution of the given
  //                                      configuration to the system's correlation
  //                                      functions.
  //-------------------------------------------------------------------------------

  template <class S, class PART, class REF>
  inline void local_correlation_functions(const int nsites,
                                          const int configs_per_site,
                                          const long int count,
                                          const int * const config,
                                          int * const work,
                                          const pcps::action_value_object<S, pcps::OpRL<S>, PART> & avoRL,
                                          pcps::CoeffValueFunc<S, PART, REF> & cvf,
                                          S * const sz,
                                          S * const sz_dot_sz,
                                          S * const s_dot_s,
                                          S * const ssss)
  {

    // check that the particles are spins
    if (typeid(PART) != typeid(pcps::Part_Spin))
      throw pcps::Exception("pcps::local_correlation_functions only implemented for spins");

    // get the spin
    const double s = 0.5 * double(configs_per_site-1);

    // initialize the work configuration
    for (int i = 0; i < nsites; i++)
      work[i] = config[i];

    // create an array to hold operators' sites' occupations
    int op_occs[2];

    // compute the contribution to the one site Sz averages
    for (int i = 0; i < nsites; i++)
      sz[i] += double(count) * pcps::unity<S>() * ( double(config[i]) - s );

    // compute the contribution to the two site Sz*Sz correlations
    for (int i = 0; i < nsites; i++)
    for (int j = 0; j < nsites; j++)
      sz_dot_sz[i*nsites+j] += double(count) * pcps::unity<S>() * ( double(config[i]) - s ) * ( double(config[j]) - s );

    // create an array to hold affected site indices
    int affected[4];

    // compute the contribution to the two site S.S correlations
    for (int i = 0; i < nsites; i++)
    for (int j = i; j < nsites; j++) {

      if (i == j) {

        s_dot_s[i*nsites+j] += double(count) * pcps::unity<S>() * s * ( s + 1.0 );

      } else {

        // initialize the value
        S v = pcps::zero<S>();

        // record which sites are changing
        affected[0] = i;
        affected[1] = j;

        // get the contribition from 0.5 * Si+Sj-
        work[i] -= 1;
        work[j] += 1;
        if ( work[i] >= 0 && work[j] < configs_per_site ) {
          op_occs[0] = work[i];
          op_occs[1] = work[j];
          v += 0.5 * cvf.get_coeff_ratio(2, affected, work, config) * avoRL.get_value(op_occs, configs_per_site);
        }

        // get the contribition from 0.5 * Sj+Si-
        work[i] += 2;
        work[j] -= 2;
        if ( work[j] >= 0 && work[i] < configs_per_site ) {
          op_occs[0] = work[j];
          op_occs[1] = work[i];
          v += 0.5 * cvf.get_coeff_ratio(2, affected, work, config) * avoRL.get_value(op_occs, configs_per_site);
        }

        // restore work configuration
        work[i] -= 1;
        work[j] += 1;

        // get the Sz contribution
        v += pcps::unity<S>() * ( double(config[i]) - s ) * ( double(config[j]) - s );

        // record the value
        s_dot_s[i*nsites+j] += double(count) * v;
        s_dot_s[j*nsites+i] += double(count) * v;

      }

    }

    // compute the contribution to the four site Si.Sj Sk.Sl correlations
    for (int i =   0; i < nsites; i++)
    for (int j = i+1; j < nsites; j++)
    for (int k =   0; k < nsites; k++)
    for (int l = k+1; l < nsites; l++) {

      // initialize the value
      S v = pcps::zero<S>();

      // get the contribitions from 0.5 * Si+Sj-SkzSlz  and  0.5 * Sj+Si-SkzSlz
      affected[0] = i;
      affected[1] = j;
      work[i] -= 1;
      work[j] += 1;
      if ( work[i] >= 0 && work[j] < configs_per_site ) {
        op_occs[0] = work[i];
        op_occs[1] = work[j];
        v += 0.5 * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(2, affected, work, config)
                 * ( double(work[k]) - s ) * ( double(work[l]) - s );
      }
      work[i] += 2;
      work[j] -= 2;
      if ( work[j] >= 0 && work[i] < configs_per_site ) {
        op_occs[0] = work[j];
        op_occs[1] = work[i];
        v += 0.5 * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(2, affected, work, config)
                 * ( double(work[k]) - s ) * ( double(work[l]) - s );
      }
      work[i] -= 1;
      work[j] += 1;

      // get the contribitions from 0.5 * SizSjzSk+Sl-  and  0.5 * SizSjzSl+Sk-
      affected[0] = k;
      affected[1] = l;
      work[k] -= 1;
      work[l] += 1;
      if ( work[k] >= 0 && work[l] < configs_per_site ) {
        op_occs[0] = work[k];
        op_occs[1] = work[l];
        v += 0.5 * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(2, affected, work, config)
                 * ( double(config[i]) - s ) * ( double(config[j]) - s );
      }
      work[k] += 2;
      work[l] -= 2;
      if ( work[l] >= 0 && work[k] < configs_per_site ) {
        op_occs[0] = work[l];
        op_occs[1] = work[k];
        v += 0.5 * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(2, affected, work, config)
                 * ( double(config[i]) - s ) * ( double(config[j]) - s );
      }
      work[k] -= 1;
      work[l] += 1;

      // get the contribitions from 0.25 * Si+Sj-Sk+Sl-
      //                            0.25 * Sj+Si-Sk+Sl-
      //                            0.25 * Si+Sj-Sl+Sk-
      //                            0.25 * Sj+Si-Sl+Sk-
      affected[0] = i;
      affected[1] = j;
      affected[2] = k;
      affected[3] = l;
      // Si+Sj-
      work[i] -= 1;
      work[j] += 1;
      if ( work[i] >= 0 && work[j] < configs_per_site ) {
        op_occs[0] = work[i];
        op_occs[1] = work[j];
        const double x = 0.25 * avoRL.get_value(op_occs, configs_per_site);
        // Sk+Sl-
        work[k] -= 1;
        work[l] += 1;
        if ( work[k] >= 0 && work[l] < configs_per_site ) {
          op_occs[0] = work[k];
          op_occs[1] = work[l];
          v += x * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(4, affected, work, config);
        }
        // Sl+Sk-
        work[k] += 2;
        work[l] -= 2;
        if ( work[l] >= 0 && work[k] < configs_per_site ) {
          op_occs[0] = work[l];
          op_occs[1] = work[k];
          v += x * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(4, affected, work, config);
        }
        work[k] -= 1;
        work[l] += 1;
      }
      // Sj+Si-
      work[i] += 2;
      work[j] -= 2;
      if ( work[j] >= 0 && work[i] < configs_per_site ) {
        op_occs[0] = work[j];
        op_occs[1] = work[i];
        const double x = 0.25 * avoRL.get_value(op_occs, configs_per_site);
        // Sk+Sl-
        work[k] -= 1;
        work[l] += 1;
        if ( work[k] >= 0 && work[l] < configs_per_site ) {
          op_occs[0] = work[k];
          op_occs[1] = work[l];
          v += x * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(4, affected, work, config);
        }
        // Sl+Sk-
        work[k] += 2;
        work[l] -= 2;
        if ( work[l] >= 0 && work[k] < configs_per_site ) {
          op_occs[0] = work[l];
          op_occs[1] = work[k];
          v += x * avoRL.get_value(op_occs, configs_per_site) * cvf.get_coeff_ratio(4, affected, work, config);
        }
        work[k] -= 1;
        work[l] += 1;
      }
      work[i] -= 1;
      work[j] += 1;

      // get the contribition from SizSjzSkzSlz
      v += pcps::unity<S>() * ( double(config[i]) - s ) * ( double(config[j]) - s ) * ( double(config[k]) - s ) * ( double(config[l]) - s );

      // record the value
      ssss[nsites*(nsites*(nsites*i+j)+k)+l] += double(count) * v;
      ssss[nsites*(nsites*(nsites*j+i)+k)+l] += double(count) * v;
      ssss[nsites*(nsites*(nsites*i+j)+l)+k] += double(count) * v;
      ssss[nsites*(nsites*(nsites*j+i)+l)+k] += double(count) * v;

    }

  }

} // end namespace pcps

#endif
